

function F = cdf_psi(wb, z, h, WB, Z)

% Returns the value of the conditional cummulative distribution function evaluated at wb given Z = z 
% Nonparametric estimation with a smoothing parameter h


%******************* Input Variables **************************************

% wb : value of the variable of interest, where the c.d.f. will be
% evaluated
% z : value of the conditional variable
% h : smoothing parameter
% WB : vector storing the outcomes of the variable of interest
% Z : vector storing the outcomes of the covariate 

%******************* Output Variable **************************************

% F : value of the cummulative distribution function

L = length(WB);
pz = size(Z,2);

% ************** kernels and auxiliary sums ******************************

numer = 0; % Auxiliary for cummulative sum
denom = 0; % Auxiliary for cummulative sum
% K = @(u) (1/sqrt(2*pi))*exp(-0.5*(u.^2));  % Gaussian
% K = @(u) ((35/32)^2)*((1-u^2)^3)*(abs(u)<=1); % Triweight
K = @(u) (3/4)*(1 - u.^2)*(abs(u)<=1); % Epanechnikov

%******** Algorithm initialisation *******************
for ell = 1:L
    
    for k=1:pz
       KF(k) = K((z(k) - Z(ell,k))/h(k+1)); % Computation of kernel function for each variable in Z
    end
    KP = prod(KF,2); % Product of kernel functions
    numer = numer + KP*(WB(ell)<=wb); % cummulative sum for numerator
    denom = denom + KP; % cummulative sum for denominator 
    
end
hz = prod(h(:,2:pz+1),2); % product of variables associated bandwidths
num=(1/(L*hz))*numer; % numerator
den=(1/(L*hz))*denom; % denominator

F = num/den; % value of cummulative distribution
end